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Summary 

Decomposition  of  a  flow  domain  into  blocks,  so-called  domain  decomposition,  is  the  most  chal¬ 
lenging  task  in  the  process  of  multi-block  grid  generation.  In  this  paper,  a  method  is  described  to 
facilitate  the  construction  of  the  blocks  in  exterior  flow  domains.  Layers  of  blocks  are  automati¬ 
cally  generated  by  solving  the  Laplace  equation  in  a  region  between  the  configuration  surface  and  a 
bounding  box  surface.  A  boundary  element  method  based  on  simple-source  distributions  is  used  to 
solve  the  Laplace  equation.  After  solving  the  Laplace  equation,  starting  with  a  face  decomposition 
of  the  configuration  surface,  blocks  are  generated  by  marching  along  the  electrostatic  field  lines 
until  a  user  specified  potential  is  reached.  The  process  may  be  repeated  to  generate  a  second  layer 
of  blocks.  However,  for  configurations  with  a  strongly  non-convex  shape,  some  blocks  may  need 
to  be  added  manually  to  prevent  deterioration  of  the  block  shapes  in  the  second-layer.  Therefore, 
the  method  is  called  semi-automatic. 
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1  Introduction 

Solving  the  Reynolds- Averaged  Navier-Stokes  (RANS)  equations  on  continuous  multi-block  struc¬ 
tured  grids  has  important  advantages  compared  to  unstructured  grids.  Concerning  efficiency,  for 
the  same  computer  memory,  multi-block  structured  flow  solvers  may  use  more  grid  cells  than  un¬ 
structured  flow  solvers  because  no  connectivity  maps  between  the  grid  cells  are  needed.  Further¬ 
more,  efficient  vectorization  and  parallelization  (on  block  level)  of  the  flow  solver  can  be  rather 
easily  obtained.  Also  multi-grid  can  be  applied  in  a  straightforward  manner  to  speed  up  the  con¬ 
vergence  rate. 

Concerning  accuracy,  a  structured  mesh  of  hexahedral  cells  is  numerically  preferable  and  especially 
suitable  within  the  boundary  layer  where  body-fitted  cells  of  large  aspect  ratio’s  (of  the  order  105) 
are  needed. 

However,  multi-block  structured  grids  are  considered  to  be  difficult  to  construct.  Especially  for 
complicated  configurations,  the  subdivision  of  a  flow  domain  in  non-overlapping  blocks  is  a  very 
labor-intensive  task  despite  the  use  of  advanced  graphical  interactive  programs. 

In  this  paper,  a  method  is  presented  to  facilitate  the  domain  decomposition  process.  The  method 
uses  concepts  as  applied  in  marching  grid  generation  techniques  like  for  example  presented  in  refer¬ 
ences  2,  3.  However,  many  of  the  problems  encountered  in  marching  (hyperbolic)  grid  generation, 
like  for  example  the  control  about  the  grid  spacing,  are  of  no  concern  because  only  blocks  are  gen¬ 
erated  and  algebraic  or  elliptic  grid  generation  methods  are  used  to  compute  the  grids  in  the  blocks 
afterwards  (Ref.  5). 

The  method  has  been  implemented  in  the  domain  decomposer  ENDOMO,  which  is  a  graphical  in¬ 
teractive  program  for  the  domain  decomposition  of  arbitrary  flow  domains  in  blocks  (Ref.  4).  EN¬ 
DOMO  is  part  of  NLR’s  Euler/Navier-Stokes  flow  simulation  system  ENFLOW  (see  Fig.  1,  Ref. 
6).  The  system  has  been  employed  to  simulate  flows  (both  on  the  basis  of  the  Euler  equations  and 
of  the  RANS  equations)  around  a  large  variety  of  complex  aerodynamic  configurations,  extending 
from  civil  and  military  aircraft  to  spacecraft.  The  grid  generator  ENGRID  is  used  to  compute  C°- 
continuous  structured  grids  for  each  block  to  obtain  final  multi-block  grids.  Given  a  C° -continuous 
grid,  the  flow  solver  ENSOLV  computes  the  solution  of  the  3D  flow  equations  based  on  Euler  or 
RANS  (Refs.  7,  8,  9). 


Fig.  1  The  NLR  multi-block  flow-simulation  system  ENFLOW 
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2  Methodology 

The  present  method  for  semi-automatic  domain  decomposition  can  probably  be  best  understood 
by  considering  the  electrostatic  field  that  would  be  generated  by  two  perfectly  conducting  closed 
surfaces  set  at  different  constant  potential  values  (Ref.  2).  One  of  the  surfaces  coincides  with  the 
configuration  surface  and  the  other  surface  is  a  bounding  box  (a  cube),  which  entirely  encloses  the 
configuration.  As  an  example,  figure  2  shows  a  half-model  of  a  generic  space-capsule  surrounded 
by  a  bounding  box.  Assume  that  the  surface  of  the  configuration  is  subdivided  by  a  non-overlapping 
set  of  faces.  Each  face  is  topologically  equivalent  with  a  unit  square  and  is  bounded  by  four  edges. 
The  points  on  the  edges  of  the  faces  can  then  be  transferred  along  the  electrostatic  field  lines  until 
a  user-defined  potential  is  reached.  In  this  way,  for  each  face  a  tube-like  block  can  be  generated. 
Such  a  block  will  then  be  bounded  by  the  face  on  the  configuration,  an  opposite  face  on  an  equipo- 
tential  surface,  and  four  faces  shaped  according  to  the  electrostatic  field  lines.  In  this  way  a  layer 
of  blocks  is  constructed  around  the  configuration  surface.  Figure  4  shows  the  in  this  way  created 
blocks  around  the  space-capsule.  This  approach  could  be  repeated  by  considering  the  outer  side 
of  the  block  layer  as  the  new  body  shape  and  by  enlarging  the  bounding  box.  Figure  5  shows  the 
enlarged  bounding  box  and  figure  6  shows  the  corresponding  far-field  block-layer.  The  method 
works  well  for  convex  body-shapes  like  a  space-capsule. 

However,  this  is  in  general  not  true  for  a  body  with  a  strongly  non-convex  shape,  like  an  aircraft,  be¬ 
cause  then  the  shape  of  the  blocks  may  deteriorate,  leading  to  unacceptable  grids  within  the  blocks. 
The  solution  to  this  problem  is  to  add  a  few  blocks  so  that  the  outer  side  of  the  block  layer  around 
the  configuration,  together  with  the  added  blocks,  define  a  more  convex  body  shape.  For  this  rea¬ 
son,  the  method  is  called  semi-automatic.  The  process  can  then  be  repeated  to  construct  a  second 
layer  of  blocks  and  also  a  third  far-field  layer.  All  blocks  together  define  a  domain  decomposition 
of  the  exterior  flow  domain. 

Figure  8  shows  a  generic  face  decomposition  of  a  wing-body  configuration.  Figure  14  shows  a 
face-decomposition  of  an  aircraft  with  T-tail.  Figure  9  and  figure  15  show  the  corresponding  first 
layer  of  blocks  around  the  configurations.  The  bounding  boxes  are  not  shown  because  again  simple 
cubes  are  used.  The  points  on  the  edges  of  the  faces  are  transferred  along  the  electrostatic  field 
lines  until  a  user-specified  height  is  reached  (instead  of  a  user-specified  potential  value  which  is  in 
general  less  suitable  for  the  creation  of  the  first  layer  of  blocks).  In  this  way  the  first  layer  of  blocks 
will  have  a  more  or  less  uniform  height.  Only  points  along  the  edges  of  the  blocks  are  stored  (wire¬ 
frame).  The  blocks  in  the  first  layer  will  contain  the  grid  boundary  layer  and  are  therefore  called 
Navier-Stokes  blocks.  The  first  layer  of  blocks  will  contain  C-type  and  O-type  blocks.  The  user- 
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specified  block  height  can  not  be  taken  to  large  because  the  shape  of  the  blocks  may  then  deteriorate 
leading  to  unacceptable  grids  within  the  blocks.  This  may  for  example  occur  near  the  junction  of 
the  trailing  edge  and  the  fuselage. 

A  few  blocks  are  added  such  that  the  outer  side  of  the  Navier-Stokes  blocks  together  with  the  added 
blocks  define  a  more  convex  body  shape.  Figure  10  shows  four  added  blocks  for  the  wing-body 
configuration.  Figure  16  shows  the  added  blocks  for  the  T-tail  configuration.  During  the  face- 
decomposition  of  the  input  configuration,  the  user  should  have  already  in  mind  that  these  blocks  can 
be  added.  This  makes  the  face-decomposition  of  a  complicated  configuration  surface  a  non-trivial 
task. 

Next,  the  process  can  be  repeated  to  generate  a  second-layer  of  blocks.  Of  course,  the  bounding 
box  must  be  enlarged.  Figure  1 1  and  figure  17  show  the  second  layer  for  the  wing-body  and  the  T- 
tail  configuration  respectively.  In  both  cases,  the  points  along  the  edges  of  the  faces  are  transferred 
along  the  electrostatic  field  lines  until  a  user-specified  potential  value  is  reached.  Thus  the  oppo¬ 
site  faces  will  lie  on  an  equipotential  surface.  The  blocks  in  the  second  layer  will  contain  no  grid 
boundary  layer  and  are  only  used  to  simulate  non-viscous  flow  features.  The  blocks  are  therefore 
so-called  Euler  blocks. 

In  general,  the  outer-side  of  the  second  layer  of  Euler  blocks  defines  a  convex  body  shape  so  that 
the  process  can  be  easily  repeated  to  generate  the  far-field  Euler  blocks.  Figure  12  and  figure  18 
show  the  far-field  blocks. 

All  blocks  together  define  a  non-overlapping  domain  decomposition  of  the  flow  domain.  Figure  13 
and  figure  19  show  the  final  domain  decompositions.  The  total  number  of  blocks  for  the  wing-body 
configuration  is  98  of  which  4  blocks  were  created  manually.  The  domain  decomposition  for  the 
T-tail  configuration  contains  168  blocks  of  which  10  blocks  were  created  manually. 

The  electric  field  can  be  computed  by  solving  the  Laplace  equation  in  the  space  between  the  body 
surface  and  the  bounding  box.  A  boundary  element  method  is  used  to  solve  the  Laplace  equation. 
Simple  source  distributions  on  the  body  surface  and  bounding  box  surface  are  used  to  represent 
the  harmonic  function.  As  an  example,  figure  3  shows  the  panel  distribution  on  the  surface  of  the 
space-capsule.  The  panel  distribution  on  the  bounding  box  surface  is  not  shown.  Details  about  the 
algorithm  are  described  in  the  following  section. 
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3  The  algorithm 

Consider  a  body  of  arbitrary  shape  and  finite  size  in  3D.  Construct  a  bounding  box  (a  cube)  around 
the  body  and  assume  that  the  interior  region  between  the  surface  of  the  body  and  the  surface  of  the 
bounding  box  is  simply  connected.  Call  this  region  V.  Next  define  the  following  potential  problem. 
Let  0  be  a  harmonic  function  satisfying  the  Laplace  equation  in  V.  Define  the  Dirichlet  boundary 
conditions  0  =  0  on  the  surface  of  the  body,  and  0  =  1  on  the  surface  of  the  bounding  box.  Then 
0  is  uniquely  defined  in  V. 

The  harmonic  function  0  is  represented  by  a  simple-layer  source  distribution  on  the  surface  of  the 
body  and  the  bounding  box,  i.e. 


=  [  k{p,q)X{q)da{q), 
Jav 


p  and  q  are  vector  variables  specifying  points  in  space  and  on  the  surface  dV  respectively  whilst 
dcr(q)  denotes  the  surface  differential  at  q  (Ref.  1).  The  kernel  function  k(p,  q)  is  defined  as 

m?)  =  I7^?ip  <2> 

and  \(q)  is  the  unknown  source  density. 

For  the  purpose  of  numerical  computation,  we  divide  dV  into  n  smooth  quadrilateral  elements 
(panels)  A j,j  —  1 ,n  and  approximate  the  function  A (q),q  €  dV,  by  a  step  function 

Ho)  =  A j,  qE  Ay,  j  =  1, . . . ,  n,  (3) 


where  the  A  j  are  some  constants. 


Correspondingly,  we  approximate  the  integral  in  Eq.  (1)  by 


m  =  ±pj  | 

J  —  J-  J 
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To  obtain  a  system  of  equations  for  the  unknown  source  densities  A j,  we  use  the  method  of  collo¬ 
cation,  applying  this  equation  at  the  midpoint  qi  of  each  quadrilateral  element  A,.  We  thus  obtain 


n 


A,-  ||  Qi~q 


;dcr(q)  =  facfr)  =  fa,  i  =  l,...,n. 


The  potential  value  fa  is  equal  0  if  A;  belongs  to  the  surface  of  the  body,  and  is  equal  1  if  At  belongs 
to  the  surface  of  the  bounding  box. 

Define  the  matrix  element  An  of  matrix  A  as 


Ali  =  Lu^7\\Mql- 


Then  we  obtain  the  system  of  linear  equations 


'y  ,  -A-ijAj  —  (f>i ,  i  —  1, . . . ,  n, 
3= 1 


or  in  matrix  notation  AA  =  <56. 


For  j  7^  i,  the  simplest  quadrature  formula  for  integration  is  used: 


^“iia-aii11^11,  i¥,i' 


where  ||  A  j  ||  is  the  area  of  A  j. 

For  j  —  i,  the  integral  in  Eq.  (6)  has  only  a  weak  (integrable)  singularity  and  can  be  integrated 
analytically  using  the  result  for  a  triangle.  If  T  is  a  triangle  of  area  ||  T  ||  and  with  side  of  lengths 
a,  b  and  c  then 


2  ||  T  ||  / a  +  b  +  c\ 

- log  -  , 

c  \a  +  b  —  cj 


(nlr) 
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where  it  is  assumed  that  pis  a  vertex  of  T,  and  c  is  the  length  of  the  edge  opposite  p  (see  Ref.  1, 
page  234).  An  is  computed  by  subdividing  A,  in  four  triangles;  each  triangle  is  formed  by  an  edge 
of  the  quadrilateral  with  q)  as  opposite  vertex.  Eq.  (9)  is  applied  for  each  triangle  and  An  is  equal 
to  the  total  sum  of  the  contributions  of  the  four  triangles. 

The  singularity  of  the  kernel  ensures  diagonal  dominance  in  the  matrix  A  and  the  problem  is  in 
general  well  conditioned  (Ref.  1).  The  solution  of  the  matrix  equation  AX  =  0  is  obtained  using 
the  subroutine  SGESV  from  the  package  LAPACK  retrieved  from  NETLIB  which  can  be  found  at 
the  Internet  address  http :  / /www.  net  lib  .  org/lapack/ index .  htral. 

The  solution  of  the  matrix  equation  gives  the  source  densities  A j,j  =  1, ...  ,n.  After  that,  the 
potential  at  a  point  p  6  V  may  be  computed  as 


and  the  gradient  of  (f>  at  p  as 


(10) 


V0(P)  =  Y,  I,  w!A^,|3  %  ~  P)  (11) 

j= l  I'  p  qi  II 

The  electric  field  is  defined  as  E  —  V</>.  Electric  field  lines  are  tangential  everywhere  to  the  electric 
field.  Starting  at  a  point  p  at  the  boundary  of  the  body,  an  electric  field  line  may  be  computed  by 
marching  in  the  direction  of  the  electric  field.  A  second-order  Runge-Kutta  procedure  is  used  to 
iteratively  calculate  the  marching  path  (Ref.  2).  A  fixed  user-specified  step  size  is  used  to  compute 
the  marching  points. 

Any  symmetry  is  of  course  taken  into  account.  Symmetry  with  respect  to  a  plane  x  —  constant,  y  = 
constant,  z  =  constant,  or  any  combinations  of  that,  is  possible.  For  example,  if  there  is  symmetry 
to  a  plane,  only  the  configuration  surface  of  a  half-model  need  to  be  panelized  and  also  only  five 
faces  of  the  bounding  box. 

For  complicated  body  shapes,  the  typical  number  of  panels  on  the  the  body  surface  and  bounding 
box  surface  is  about  10000  and  1000  respectively,  so  that  the  matrix  A  contains  about  108  elements. 
The  matrix  equation  is  solved  on  a  SX5  supercomputer.  The  corresponding  CPU  time  is  about  5-10 
minutes. 


4  Conclusions 


A  method  has  been  developed  to  compute  block  layers  about  arbitrary  configurations.  For  strongly 
non-convex  configurations,  additional  blocks  must  be  added  to  prevent  deterioration  of  the  block 
shapes.  After  that,  second  or  third  layers  of  blocks  are  generated  automatically,  leading  to  a  satis¬ 
factory  domain  decomposition  of  the  whole  flow  domain.  The  method  appears  to  be  an  important 
tool  to  reduce  the  human-time  to  generate  multi-block  grids  for  general  3D  configurations.  The 
user-time  to  complete  the  domain-decomposition  for  a  simple  geometry,  like  a  space-capsule,  takes 
a  few  hours.  For  a  complicated  geometry,  like  an  aircraft  with  T-tail,  it  takes  a  few  working  days.  In 
general,  the  time  needed  for  domain  decomposition  has  become  comparable  with  the  time  needed 
for  surface  modelling  (i.e.  making  a  CAD  geometry  suitable  for  CFD  computations). 
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Fig.  8  Generic  face-decomposition  of  wing-  Fig.  9  First  layer  of  automatically  created 


body  configuration 


blocks 


Fig.  10  Four  blocks  added  to  make  outer 


shape  more  convex 


Fig.  12  Third  layer  of  automatically  created 


Fig.  1 1  Second  layer  of  automatically  created 


blocks 


Fig.  13  Final  domain  decomposition 


far-field  blocks 
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Fig.  14  Face  decomposition  of  an  aircraft  with  Fig.  15  First  layer  of  automatically  created 
T-tail  blocks 


Fig.  16  Blocks  added  to  make  outer  shape  Fig.  17  Second  layer  of  automatically  created 
more  convex  blocks 


Fig.  18  Third  layer  of  automatically  created  Fig.  19  Final  domain  decomposition 
far-field  blocks 


